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ABSTRACT 



This thesis examines two implementations of the para- 
bolic equation approximation to the acoustic wave equation 
aimed at removing three errors inherent to the wide-angle 
parabolic equation (WAPE) model. First, the selection of 
the range-step size used by the split-step Fourier algorithm 
affects the convergence of the solution. Second, in certain 
ocean environments WAPE incorrectly computes the down-range 
transmission loss. Finally, WAPE does not reproduce the 
standard normal mode basis set as defined by normal mode the- 
ory. A double-precision implementation of the WAPE (DP- 
WAPE) is developed to evaluate the dependence of solution 
convergence on the numerical precision of the model. 

Finally, an implementation that is insensitive to the choice 
of the reference sound speed (COIPE) is evaluated for its 
ability to reduce or remove the latter two of these three 
errors. The stability of the WAPE solution was found to be 
unaffected by the DP-WAPE implementation. The range-step 
dependence is inherent to the split- step algorithm. The 
COIPE corrects the transmission loss anomaly and satisfacto- 
rily reproduces the standard normal mode basis set. 
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I. INTRODUCTION 



The fall of the former Soviet Union and the resulting 
shift in the international political power base has resulted 
in a dramatic shift in both the content and the application 
of United States Naval forces during the past decade. The 
world's largest blue water navy finds itself increasingly 
involved in green water operations. This shift to littoral 
waters has brought to the forefront the ever present danger 
of easily obtainable and extremely capable diesel submarines 
and shallow water mines. Countries with otherwise limited 
military resources can acquire the ability to thwart or at 
the least hinder the application of national objectives. As 
a result of these changes, numerous efforts toward the 
development of more accurate acoustic propagation prediction 
models have arisen. These include areas such as matched- 
field processing, transient localization, and underwater 
acoustic communications. Essential to these efforts is the 
development of acoustic prediction models which can produce 
both valid and accurate results in highly variable shallow 
water ocean environments . 
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A. OCEAN ACOUSTIC MODELING TECHNIQUES 

Current acoustic modeling efforts generally fall into 
one of four broad categories. These are ray models, wave- 
number integration techniques, normal mode models, and para- 
bolic equation models. Each of these are based on a variety 
of approximations to the linear acoustic wave equation. 

Many of the early models are based on the geometrical limit 
and define ray trajectories of sound propagation. The 
resulting ray models can quickly predict pulse propagation 
travel times and provide rough approximations of the sound 
field. The limitation of this approach is its reliance on a 
high-frequency limit which neglects finite- frequency effects 
such as diffraction. In addition, it has been shown (Smith, 
et. al., 1992) that the ray equations form a set of coupled, 
nonlinear equations which suffer from chaotic solutions in 
the presence of range-dependence. Given the current impor- 
tance of low-frequency propagation in range-dependent envi- 
ronments, the usefulness of ray models in shallow water 
media is limited. 

Wavenumber integration methods produce highly accurate, 
full-wave solutions to the wave equation. However, they are 
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computationaly intensive and are, by design, limited to 
range- independent environments. They can be quite useful in 
determining near-field effects but are not of practical 
utility in predictions of far-field acoustic propagation. 

Normal mode models are based on a separation of vari- 
ables in range and depth of the acoustic wave equation. 

This approach maintains the full-wave characteristics of the 
acoustic field, and provides highly accurate results to both 
simple and complex range dependent ocean environments . The 
solutions produced are based on the single-frequency Helm- 
holtz wave equation. Due to the computational complexity 
involved with high frequencies and broadband pulses, as well 
as range-dependent media, such models are best suited for 
low-frequency, shallow water environments which support only 
a small number of propagating modes. The results remain very 
accurate, but computational time can become counter produc- 
tive . 

The final acoustic modeling technique, and undoubtedly 
the most popular in range-dependent environments, is based 
on the parabolic approximation to the Helmholtz wave equa- 
tion. Numerous marching algorithms have been developed that 
solve the acoustic field as a boundary value problem. While 
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it is a one-way wave equation treatment, it is a full-wave 



model and therefore includes all of the finite- frequency 
effects of diffraction. Two-way adaptations of parabolic- 
equation models have been developed, but do not result in 
significant changes for the majority of shallow water envi- 
ronments of interest. Solutions are, however, still single 
frequency, requiring multiple solution generation for pulse 
propagation responses. It is this method that will be the 
primary focus of this thesis. 

B. HISTORY OF THE PARABOLIC EQUATION APPROXIMATION 

The parabolic equation (PE) approximation method was 
first applied to the problem of underwater acoustics by Tap- 
pert (1977) . Its historical roots are much deeper, however. 
Breakthroughs in physics and mathematics in the mid-1920's 
provided the basis for the parabolic wave equation used in 
acoustics today. The standard parabolic equation (SPE) has 
the same form as the Schrodinger ' s equation in quantum 
mechanics. The earliest documented use of the SPE as an 
approximation to the theory of wave propagation dates to the 
mid-1940's and the work of Leontovich and Fock (1946) who 
originally coined the term "parabolic equation method" . 
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They applied the method to predicting the diffraction on 
long range tropospheric wave propagation caused by the 
spherical shape of the earth. Fock (1965) expanded this 
approach to include high-frequency scattering as well as 
microwave propagation in waveguides . 

With the advent of lasers and their coherent radiation 
sources in the 1960 's, the PE method was further extended to 
laser beam propagation where it is generally referred to as 
the "quasi-optical" equation. The quasi-optical equation is 
most widely used in nonlinear optics where the index of 
refraction depends on intensity. It is often called the non- 
linear Schrodinger ' s equation. The PE method has also been 
applied, in more recent years, to the field of fiber optics 
and plasma physics . 

Beam propagation in random media also lends itself well 
to the PE method. Regardless of the source of the beam -- 
radio, acoustic, optical -- the problem is analogous to the 
quantum mechanics problem of motion in a random potential . 
The radar application of this problem in a randomly fluctu- 
ating ionosphere using the split-step Fourier algorithm, 
discussed later in the next chapter, is covered extensively 
by Hardin and Tapper t (1974) . 
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The most extensive use of the PE approximation to the 
elliptical wave equation has been in the field of low-fre- 
quency underwater acoustic propagation. More than 120 arti- 
cles and technical reports regarding the PE method and its 
applications to underwater acoustics were published in the 
last 15 years. Numerous marching algorithms and computer 
applications have been constructed with varying degrees of 
success. The earliest of these were based on the split-step 
Fourier algorithm (Hardin and Tappert, 1974) . In general, 
these models accept input defining the sound speed, volume 
loss profiles, and depth contours, and produce output of the 
acoustic field and associated transmission loss (TL) curves. 
These results are generally in excellent agreement with 
experimental measurements for deep ocean environments. The 
theory behind this method and subsequent attempts to improve 
it are covered in detail in the next chapter. Since its ini- 
tial use, the PE method has been extended to include higher 
acoustic frequencies, random internal-wave fluctuations in 
the index of refraction, and shallow-water environments with 
mixed levels of success. 
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C. ADVANTAGES AND LIMITATIONS OF THE PE METHOD 

Long-range, low-frequency sound propagation is gener- 
ally dominated by rays having small grazing angles since 
rays propagating at steep angles are greatly attenuated due 
to penetration and absorption in the seabed. Tappert intro- 
duced the PE method, which decomposes an elliptical wave 
equation into two equations through the separation of outgo- 
ing and incoming propagating fields. Tappert 's implementa- 
tion of this method was the first applied to underwater 
acoustics and is referred to as the standard parabolic equa- 
tion (SPE) . The SPE implementation and many of its more 
recent incarnations share a number of inherent advantages 
and limitations. 

The advantages of the PE approximation can be divided 
into three categories (McDaniel, 1975) . First, in the long- 
range, low-frequency environment outlined above, the govern- 
ing equation can be easily expressed by a parabolic equa- 
tion, which is easier to solve than either the elliptical or 
the hyperbolic types. The PE method, therefore, provides a 
relatively quick solution to long-range, low-frequency, 
range- dependent problems. Second, in solving the Helmholtz 
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equation, the problem is posed as a pure boundary value prob- 
lem. In practice, the determination of the vertical bound- 
ary condition is often difficult and imprecise. The PE, 
however, poses the problem as one of an initial -boundary 
problem, therefore avoiding the difficulty of the vertical 
far-end wall boundary condition. Finally, the PE solution 
is generally cheaper to obtain in both memory and run-time 
than solving the elliptical equation. 

The limitations of the PE method can be placed in one of 
two broad types. The first of these results from the mathe- 
matical formulation of the PE itself which does not allow the 
inclusion of certain physical phenomenon. The PE is only 
valid in the far- field and therefore cannot be used for a 
complete analysis of near-field effects. The index of 
refraction is assumed to be slowly varying in range. Large 
and or abrupt changes in the index of refraction can limit 
the accuracy and validity of results. The PE uses the one- 
way wave equation and wave propagation from the reverse 
direction is ignored. Backscattering effects are, there- 
fore, neglected. Treatment of the so-called "square root 
operator" defines the size of the propagation angle which in 
turn effects the validity of the PE. 
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The second type of limitation is defined by the method 
of solution to the PE. Many of the limitations of this type 
are specific to the solution algorithm being implemented. 

The fast Fourier transform works well only for equations 
with constant coefficients, and therefore must be used with 
caution. If the boundary is rigid, the algorithm cannot 
treat it exactly. Finite difference techniques require reg- 
ular grid partitions. Using the horizontal-interface PE to 
handle irregular interfaces, unless the sloping angle is 
small or density change is small, may not produce satisfac- 
tory results. In theory, there are no limitations on fre- 
quency. In practice, however, using finite difference 
techniques for high-frequency problems results in intolera- 
ble and expensive computational times. 

The two types of limitations result in three basic 
errors. The first error is due to the initial PE approxima- 
tion. Findings by Smith and Smith (1997) and others show 
that in the SPE a normal mode will be propagated with the 
correct amplitude and mode shape but with an error in phase 
velocity. Some higher order PE approximations improve the 
accuracy of the phase velocity, but may not properly propa- 
gate the mode amplitudes. The error in phase velocity is 
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addressed in a number of PE approximation methods in the way 
in which the square root operator is expressed. The second 
error is introduced by the selection of the range step size. 
The error is typically found in the second or third order of 
the range step increment (Jensen, et. al., 1994). If the 
error is defined as E n/ then for the case of several modes 
propagating, the error produced in a single range step size 
would be E = ^E n . The third type of error is due to trun- 

n 

eating the field. The initial approximation of the PE poses 
basic limitations as outlined above. The errors introduced 
by truncating the field are rooted mainly in computational 
efficiency in which the numerical results are obtained. A 
detailed analysis of the effects of range step choice were 
investigated by Smith (2000) . 

D. THESIS SUMMARY 

The objective of this thesis is to explore the accuracy 
and validity of various approaches to the PE approximation 
method as it applies to underwater acoustics, in particular, 
how each of these methods or implementations affects the 
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errors outlined above. The resulting accuracy and validity 
of the model output will be analyzed. 

Chapter II, Theoretical Background, reviews the theory 
relevant to the material in the subsequent chapters . Over- 
views are provided of the parabolic equation approximation 
used for acoustic propagation modeling, the standard para- 
bolic equation (SPE) method, the wide angle parabolic equa- 
tion method (WAPE) , and the reference sound speed (c Q ) 
insensitive parabolic equation method (COIPE) implementa- 
tions. Additionally, a brief overview of the method of nor- 
mal modes is discussed with emphasis on the modal 
decomposition of acoustic fields generated by PE models. 

Chapter III, Numerical Implementation, outlines the 
numerical modeling techniques used in implementing the c Q - 
insensitive and double-precision WAPE models. 

Chapter IV, Numerical Results, explores the results 
obtained from each of the implementations examined. These 
include SPE, WAPE (single and double precision) , and COIPE. 
Each model is examined for validity and accuracy against a 
number of standard test cases. The models were also analyzed 
for their ability to reproduce the standard normal mode 
basis set. 
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Chapter V, Conclusions and Summary, presents a 



summa- 



tion of conclusions and identifies areas of further investi- 
gation and research. 
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II. THEORETICAL BACKGROUND 



This chapter outlines the theoretical background sup- 
porting the development of the standard, wide-angle, and c 0 - 
insensitive implementations of the parabolic equation (PE) 
approximation to the acoustic wave equation. The bulk of 
this theoretical development is directed toward its applica- 
tion in the Mont erey-Mi ami PE (MMPE) implementation of the 
wide-angle parabolic approximation. The MMPE implementation 
is emphasized here since it has been subjected to extensive 
analysis and numerous numerical improvements over its life 
span. It also serves as the framework of the COIPE model 
discussed and analyzed later. A summary of the method of 
normal modes with an emphasis on the modal decomposition of 
acoustic fields generated by PE models is also included. 

The entry point in deriving the PE is the Helmholtz 
equation, in cylindrical coordinates. 



where p(r,z) is the acoustic pressure, k 0 = 0)/c 0 is the ref- 
erence wavenumber, n(r, z, (p) = c 0 /c(r, z, (p) is the acoustic index 
of refraction, c 0 is the reference sound speed, and c(r, z, (p) 




( 2 . 1 ) 
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is the acoustic sound speed. All of the environmental char- 
acteristics are represented in n(r, z, cp) . The treatment of 
density contrasts will not be presented here, but is 
included in the MMPE model by an effective index of refrac- 
tion term (Smith, 2000) . A time -harmonic component of the 
acoustic field then has the form 

P(r, z, (p, cot) = p(r, z, cp)e _10)t . (2.2) 

A. OPERATOR NOTATION 



In order to develop the various parabolic approxima- 
tions to the Helmholtz equation, it is useful to introduce an 
operator notation. Let 



P 



op 



a_ 

9r 



(2.3) 



and 



Q op = (H + e + u+l) 1/2 , 



(2.4) 



, 2. 1 a 2 . i a 2 „ , 

where E = n z - 1 , jx = — , and x> = * To account f or 

the dominate cylindrical spreading and to simplify the form 
of the Helmholtz equation, the acoustic pressure may be 



defined as p(r, z) = - 7 q(r, z) , which yields 
Jr 
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( P op + k oQo p )q = 0- (2.5) 

Proper factorization of this equation is accomplished by 



defining q = Q“jJ /2 u (Tappert, 1977) , such that 



(Pop + ik 0QopK P op- ik 0Qop) u + ik 0[ P op>QopJ u = °- 



( 2 . 6 ) 



The commutator [P op ,Q op ] is exactly zero for layered media, 
and is assumed negligible. The remaining factors of Eq. 
(2.6) define the wave equations for the incoming and outgo- 
ing fields, the latter of which then satisfies 



In cases where backscattered energy may be neglected, Eq. 
(2.8) provides the full description of the forward propagat- 
ing acoustic energy. 

B. THE SPLIT-STEP FOURIER ALGORITHM 

The most common methods of computing PE solutions are 
(1) the finite element method (Collins, 1994), (2) the 

implicit finite difference method (Lee, et. al . , 1981), and 
(3) the split-step Fourier (PE/SSF) method (Hardin and Tap- 



P op u = ik oQop u 



(2.7) 



or 




( 2 . 8 ) 
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pert, 1973) . Of these three methods, the PE/SSF method has 
the distinct advantage of speed. This is accomplished by- 
decomposing the acoustic field into a slowly modulating 
envelope function and a phase term which oscillates with the 

acoustic frequency. Defining u = and substituting into 

Eq. (2.8) yields 

= -ik 0 ¥ + ik 0 Q op 'F = -ik 0 H op 'P , (2.9) 



where H op = l-Q op is a Hamiltonian-like operator which 
defines the evolution of the PE field function in range. 

In order to compute the solution to the acoustic field, 
a marching algorithm must be developed of the form 

^(r + Ar) = $>(!■)¥( r) . (2.10) 

The propagator, 3>(r) , is designed to propagate the solution 
out in range. The PE/SSF method accomplishes this via 

O(r) = e -ik ° Hop(r ) Ar (2.11) 

where 



Hop(r) 



i f (r+Ar) 

ArJ H op( r ') dr • 

r 



( 2 . 12 ) 



In practice, the Hamiltonian is usually considered constant 



over a small range-step such that H op (r) = H op (r) . 
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The operator H op (and Q op ) is not a scalar operator but 



a combination of scalar and differential operators. Each 
individual operator within H op can, however, be applied by a 
simple scalar multiplication in the proper domain. There- 
fore, each of the terms of H op must be separated in order to 
comply with the PE/SSF algorithm. This necessitates the 
approximation of the square-root operator, specifically 

HjX™n(r > z)VTj^ + U 0D (n(r f z)) + vJ^) . (2.13) 



op Wz 2 * r 2 3(|) 2 ’ 



op Uz 2 



OI \r 2 a<t> 2 J 



It is this approximation that forms the basis of numerous 
higher-order PE model implementations. 

If the uncoupled azimuth approximation is employed 
then V op = 0 identically. The U op is a multiplication opera- 
tor in z-space and is therefore a diagonal matrix. The oper- 
ator T op is not diagonal in z-space, resulting in different 
eigenfunctions being coupled. The T op is, however, diagonal 
in vertical wavenumber space. It is advantageous then to 
treat each operator separately, one in k z -space, and the 
other in z-space. Employing the Baker-Campbell-Hausdorf f 
expansion (Bellman, 1964), 

e A + B _ e A e B e [A,B]+[A,[A,B]] + (B,(B,A]] + ... # (2.14) 
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T op and U op are both 



where A = -ik 0 ArT op and B = -ik 0 ArU op . 
small, therefore it is assumed their products of second 
order are negligible. Using Eq. (2.14), <I>(r) becomes 



. _ik oT U op( r + Ar ) -ik fl ArT ~ ik oT U op( r ) _ 

O(r) = e 2 e 1K o Ar, o Pe 2 (2.15) 

Third order accuracy results from this centered step scheme 
(Jensen, et. al., 1994). Due to the formulation of the prop- 
agator there are no intrinsic losses due to the numerical 
scheme and energy conservation is retained (Tappert, 1995) . 

In summary, the PE/SSF algorithm accepts the PE field 
function *¥ specified at some range r in the z-domain. This 



is followed by a multiplication of e 2 , the z-space 

operator defined at the beginning of the range-step. This 
product is then transformed to the k z -domain and multiplied 



by e lk ° ArT °p( k z) , the k z -space operator. Another transformation 
to the z-domain and multiplication by the z-space operator 

-ik 0 ^U op (r + Ar) 

e 2 , defined at the end of the range-step, produces 

the final resultant field function at r+Ar in the z-domain. 
Note that the FFT convention used here and in the numerical 
implementations is 
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'i'(z) = FFT(^(k z )) # (2.16) 

*F(k 7 ) = IFFT( x F(z)) ‘ 

Wrapping these steps into a single equation results in 



x F(r + Ar,z) = 



(2.17) 



-ik 0 yU op (r + Ar,z) 

e 1 



x FFT< e - ik o ArT °p( k z) x 




e -'k 0 yU op (r 1 z) x 



^(r, Z)] 1 • 



C. STANDARD PE DEVELOPMENT 

The lowest order approximation to the Hamiltonian oper- 
ator is commonly referred to as the standard parabolic equa- 
tion (SPE) . This approximation is obtained by assuming the 
operators are small, 



The first of these is merely a recognition that sound speeds 
vary by less than 2% in most ocean environments, so n(r, z) = 1 
everywhere. The second assumption can be seen to imply that 



accuracy to small angles. 

With these approximations, the SPE is obtained by a 
binomial expansion of the square root operator. 



e « 1 and « 1 



(2.18) 



vertical angles = sin 2 G are small. This limits the SPE 

l k 0 J 



fk? 
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(2.19) 



H, 



spe 




The separated operators then have the form 




( 2 . 20 ) 



and 



U s P e( r ’ z ) = -|(n 2 (r,z)- 1) . 



( 2 . 21 ) 



Note that in the context of H spe acting as a Hamiltonian 
operator, the operators T spe and U spe correspond to kinetic 
and potential energy operators, respectively. In the verti- 
cal wavenumber domain. 



D. WIDE-ANGLE PE DEVELOPMENT 

A higher order approximation to the operators was 
introduced by Thompson and Chapman (1983) and is often 
referred to as the wide-angle PE (WAPE) approximation. The 
operators of the WAPE are defined by 




( 2 . 22 ) 




(2.23) 
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and 



U wape( r ’ z ) = -(n(r,z)-l) . 



( 2 . 24 ) 



In wavenumber space, the WAPE kinetic energy operator takes 
the form of 



Note that energy propagating at angles beyond vertical, 
k z >k 0 , is evanescent since 



The wide-angle PE, and other higher-order approxima- 
tions to the SPE, was originally developed in an effort to 
extend the region of validity of the SPE. While the wide- 
angle PE succeeded in this goal, it introduced a number of 
other issues. These include a sensitivity to the selection 
of the reference sound speed and an inability to reproduce 
the standard normal mode basis set. This said, however, the 
wide-angle approximation has been very successful in produc- 
ing useful and valid results in many situations. The follow- 
ing section addresses an attempt to mitigate the reference 
sound speed sensitivity experienced by the WAPE approxima- 
tion. 




( 2 . 25 ) 




( 2 . 26 ) 
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E . COIPE DEVELOPMENT 



WAPE and other higher-order models are not exact repre- 
sentations and may under certain circumstances produce 
results worse than the SPE which they are meant to improve 
upon (Porter and Jensen, 1993) . In these cases, it is often 
found that the error results from errors in the underlying 
normal mode basis set, their associated phase speeds, and a 
sensitivity to the choice of reference sound speed, c 0 . The 
choice of c 0 is the one ambiguous feature of all PE models. 
Due to this ambiguity, it is desirable to develop a model 
that is insensitive to the input reference sound speed. 

Then, any significant deviations in the choice of c 0 would 
not adversely affect the computed results. In an attempt to 
develop a wide angle implementation that is insensitive to 
c 0 , Tappert (1991) developed the c 0 -insensitive PE (COIPE) 
approximation . 

The COIPE approximation based on the WAPE approximation 
with the operator replacement 



1 3 1 3 ,, N 1 3 

‘ — n-‘(r,z)- 



kga z 2 k o (r)0z v ’ ; k 0 (r)3z ' 



(2.27) 
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where k 0 (r) is now range -dependent due to the range- dependent 



nature of the reference sound speed, c 0 (r) . Introducing the 
operator 



the c 0 -insensitive approximation to the Hamiltonian operator 
becomes 



In this form the spherical wavefront property of the TC 
approximation is maintained for steep angles. Furthermore, 
this approximation is not sensitive to the choice of c 0/ and 
reduces to a c 0 - independent approximation for small angles 
(Tappert, 1977) . The resulting wave equation approximation 
is 



To implement this approach, a transform to a new depth 
coordinate is necessary. The transformation relationships 
are 




(2.28) 



H Ins (r,z) = ( 1 - n(r, z)) + ( 1 - J~\ - p(r)n-> (r, z)p(r)) . 



(2.29) 



57 = -ik 0 H lns (r)'F . 



(2.30) 



Z 




(2.31) 



0 
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n(z, r) = n(z, r) , 



(2.32) 



and 

'F(z,r) = n- I/4 (z,r)'I'(z > r) . (2.33) 

The corresponding Hamiltonian operator then becomes 

H(r) = (l-n(z,r)) + (l--7l-p 2 (r)) . (2.34) 

Note that in the tilde-domain this has the same form as the 
TC operator. The resulting form of the parabolic wave equa- 
tion in the tilde domain is 

p¥ = ik 0 Hms(r)'F . (2.35) 

or 

The corresponding kinetic and potential energy operators of 
the COIPE in the tilde-domain may then be represented by 

and 

Uins(r, z) = 1 - n(r, z) . (2.37) 

The COIPE requires, in general, the computation of c 0 
for each range step. To compute c 0 , the condition that the 
water depth be unchanged in the z coordinate must be satis- 
fied. This is accomplished via 
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(2.38) 



h(r) 



h(r) = J */n(z, r)dz , 



o 



where h(r) is the water depth at range r. Under this condi- 
tion, the reference sound speed is 



The University of Miami PE/SSF acoustic model (UMPE) 
has been tested using the technique outlined above with 
excellent results (Tappert, 1995) . A numerical implementa- 
tion of the COIPE within the existing framework of the MMPE 
model is covered in the following chapter. Comparative 
analysis of this approach with the SPE, MMPE, and selected 
normal mode test cases is covered in Chapter IV. 

F. NORMAL MODE FUNCTIONS AND THEIR DECOMPOSITION FOR THE 
PE APPROXIMATION 

One of the difficulties encountered with the WAPE 
approximation is that it does not decompose into the stan- 
dard normal mode basis set as defined by normal mode theory 
and the SPE. Thomson (1993) looked at this problem for the 
SPE approximation. He showed that the field p that satisfies 




0 



(2.39) 
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the acoustic wave equation exactly relates to the field *F 
satisfying the SPE approximation. The wavenumber of 'F and 
the modal amplitudes may be determined exactly. The corre- 
sponding normal mode amplitudes and wavenumbers are then 
obtained. The WAPE, however, has modal eigenfunctions and 
eigenvalues of 'F distinct from the acoustic wave equation. 

Smith and Smith (1997) showed that the modal decomposi- 
tion of solutions based on the WAPE generated erroneously 
range-dependent mode amplitudes in a range- independent envi- 
ronment. This was due to the mismatch between the standard 
normal mode basis set (based on the Helmholtz equation) and 
the proper basis set of the WAPE. Furthermore, they were 
able to develop an approximate numerical scheme for comput- 
ing the proper WAPE basis set. This resulted in the proper 
range- independent character of the mode amplitudes. Addi- 
tionally, they suggested that the higher-order COIPE version 
of the WAPE might overcome such issues and decompose prop- 
erly into the standard normal mode basis set. This will be 
one of the primary tests of the improved accuracy of the 
COIPE over the WAPE. 

Normal modes are obtained by using the method of sepa- 



ration of variables. As presented here, it assumes the ocean 
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is horizontally stratified, is of constant water density, is 
azimuthally symmetric, and is range-independent. As was 



presented for the PE, the time harmonic acoustic field 
p(r, z)e -ia)t at r>0 due to a point source at z = z s and r = 0 sat- 
isfies the 2D Helmholtz equation (no azimuthal variations 
possible) 



is obtained. Here, , the square of the horizontal wave- 
number for mode m, is the separation constant (eigenvalue) 
and v m (z) is the specific mode function (eigenfunction) 
associated with the separation constant. Given an arbitrary 
pressure field, the sum of the normal modes is required and 
is expressed as 



2 




(2.40) 



For solutions of the form p(r, z) = -y<|)(r)v(z) , the depth- 

Vr 



dependent modal equation 




(2.41) 



P( r - Z > = -J T £ 4 > m( r ) v m( z ) ' 



(2.42) 



m = 1 
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where the cylindrical spreading factor has been explicitly 



included. The factor <j> m (r) can easily be shown to have the 

form of a Hankel function scaled by mode excitation values 
depending on the source strength and depth in the waveguide. 
In the far-field, the only range-dependence of this function 
is in the phase (having already accounted for cylindrical 
spreading). Thus, a generalized mode amplitude, A m (r) , may 
be defined according to 



P( r ,z) = I^A m (r)v m (z) . 



( 2 . 43 ) 



The normalized eigenfunctions form a complete, orthogo- 
nal basis set, defined by 
D 



1 

p(z) 



Jv„(z)v m (z)dz = 5, 
0 



( 2 . 44 ) 



Thus, by using the orthogonality relation and assuming 
normalized modes, the generalized modal amplitudes can be 
extracted from the computed pressure field by forming an 
inner-product 

D 

. ,, r fp( r ’ z ) v m( z )j _ 

A "> (r) = /r J p(z) ' dZ - (2 ’ 45 

0 
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This technique is used in the analysis and comparisons of the 
normal mode amplitudes of the SPE ; WAPE, and COIPE models in 
Chapter IV. 

Finally, the corresponding depth separated equations of 
the various parabolic approximations introduced should be 
examined. First, note that the normal mode equation, Eq. 
(2.41), may be written in terms of the operators |J, and £ 
(defined in Eq. (2.4)) as 



This indicates that the basis set, v m , are eigenfunctions of 

the combined operator (n + e) with eigenfunctions (K 2 m /k 2 0 -l). 

The PE approximation has a depth separated equation of 
the form 



where the functions % m are the eigenfunctions of the corre- 
sponding PE approximation. In order for the eigenfunctions 
X m and v m to be the same (up to some normalization factors) , 
the operator H op must be some rational function of the oper- 
ator (n+e). For example, a simple polynomial series with 




(2.46) 




(2.47) 
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terms (ji + s) n , where n is some integer, will have the same 



basis set v m . 

The standard parabolic equation is a simple first order 
binomial expansion of the Hamiltonian operator, i.e. 

H sp e = -^ + 8) ' (2 ' 48) 

This obviously satisfies the above criterion, therefore the 
SPE approximation has the same normal mode basis set as the 
fundamental Helmholtz equation. 

The wide-angle PE, on the other hand, may be written in 
the form 

H wape = 2-VTTjI- 7T+E • (2.49) 

Performing a power series expansion on this operator to sec- 
ond order yields 

H wape=— ^ + 8 ) + §(H 2 + e 2 ) • (2.50) 

This shows that the WAPE approximation is only first order 
accurate in the operator (^i + s) but has errors in the second 
order terms. It is this weakness of the WAPE which generates 
some of the associated problems with this approximation. 

The c 0 - insensitive PE approximation, however, may be 
written as 
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(2.51) 



H 



Ins 



= 2- /i + -L-VTTi. 

«/ jm 



A power series expansion of this operator to third order 
yields 



H ins 88 - + e) + + e) 2 - 17(H- + e ) 3 + ^ • 



16 



16 



(2.52) 



Thus, the COIPE is accurate to second order in (fi + s) with 
leading error in the third order terms. This was first shown 
by Tappert and Brown (1996) . It is this feature of the COIPE 
which makes it an attractive alternative to the WAPE approx- 
imation. As will be shown, this also generates solutions 
which may be decomposed into the standard normal mode basis 
set without significant error. 
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Ill . NUMERICAL IMPLEMENTATION 



This chapter outlines the numerical implementation of 
the PE models used in the analysis of the c 0 -insensitive PE 
approximation (COIPE) as well as a double-precision version 
of the MMPE. The MMPE (Smith, 2000) is an existing well-doc- 
umented version of the WAPE approximation. In the cases 
where a SPE solution was required for comparison purposes, 

the operators U wape and Twape were replaced with U spe and Tspe 
within the framework of the MMPE model. As was discussed in 
the previous chapter, the COIPE implementation strives to 
eliminate the sensitivity to the chosen reference sound 
speed suffered by the WAPE approximation. As will be shown 
in the next chapter, a small change in the specified refer- 
ence sound speed can result in a dramatic change in the pre- 
dicted acoustic field. This error often increases with 
range and for complex, highly range -dependent environments 
can quickly result in the produced data being unusable. The 
numerical implementation described here was chosen such that 
it could be easily integrated into the existing MMPE model to 
provide the most reliable comparison of results. The 
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results of these comparisons are discussed in the following 



chapter . 

A. ENVIRONMENTAL GRID SIZES 

The MMPE , as well as the COIPE, implementation dis- 
cussed here uses a discretized grid to describe the environ- 
ment. The grid step-sizes are user definable and can 
dramatically affect the solution obtained by the model. The 
grid sizing method described here is that used by the MMPE 
model as well as the COIPE implementation of the MMPE. The 
environment is discretized by a mesh of size (Ar, Az) . This 
results in the field and propagator functions being dis- 
cretized in depth with array length N, where N is the corre- 
sponding length of the FFT used in the PE/SSF algorithm. In 
other words. 



'F(r,z)=>'F(r i , z n ) 



( 3 . 1 ) 



where 




( 3 . 2 ) 



and 
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The depth mesh is defined such that grid points lie on frac- 
tional values of Az . This convention results in the avoid- 
ance of carrying the zero-pressure value at z = 0 throughout 
the calculations. It should be noted that because of using 
the full Fourier transform, half of the depth mesh values 
define an image ocean for negative depths. This has the 
added benefit of enforcing the surface boundary condition, 
defined in the next section, via FFT symmetry. Analysis by 
Tappert (1977) and Smith (2000) showed that by considering a 
lower limit on the allowable angles of propagation, a 
default value of Az may be defined and therefore a default 
transform size N. This is possible since the depth mesh 
influences the wavenumber increments Ak z via the FFT. For a 

given environment with a relatively small angle of propaga- 
tion, the mesh size may be increased. An increase in the 
mesh size results in a corresponding decrease in the compu- 
tational time. As propagation angles increase, the mesh 
size is reduced. The most accurate solution, therefore, 
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should be obtained when both Az and Ar are on the order of a 
wavelength. 

B. TILDE TRANSFORMATION AND C 0 GENERATION 



The tilde transformation begins by defining a rescaled 
depth variable z in terms of the true depth z 



JVn(z', r)dz' . 



( 3 . 4 ) 



If it is required that the bottom depth remains fixed, in 
other words that Zb = z b = h , then the range -dependent bottom 
depth is 

h(r) h(r) h(r) 



h(r) = J Vn(z\ r)dz’ = J J-j^-dz’ = J^T) f —== 

J J V c(z,r) ^ u J Jc(z\ 1 



:dz’ 



( 3 . 5 ) 



0 0 0 
where c 0 has been specifically declared as range-dependent, 
due either to h(r) or c(z,r). Solving for c 0 (r), 



c 0 ( r ) = 



h(r) 



-i-2 



-L f_L. 

h( r ) J Jc(z\ r) 



dz’ 



( 3 . 6 ) 



This expression defines the rang- dependent reference sound 
speed to be used in the parabolic approximation. Its decla- 
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ration then eliminates the need for a user-defined value. 
Furthermore, since the definition is based on an evaluation 
of the local environment at each range, it will preserve rec- 
iprocity. 

If c(z) varies linearly with depth. 



c(z) = c(zj_ j) + b(z -Zj_ j) = a + bz. 



(3.7) 



where b is the sound speed gradient, 



c(z) — C(Z; _ i ) 



(3.8) 



and 



a = c ( z j _ i ) — bzj _ , , 



(3.9) 



then the integral may be written 




(3.10) 



Using a table of integrals, this becomes 




(3.11) 



Substituting into Eq. (3.11) for a and b. 




(3.12) 
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Now consider the sampling of the environment within the 
model. The typical definitions are 



z k = ( k -|) Az ' 

k = 1,2,3...* ' 



(3.13) 



where N is the FFT size. Based on this grid sampling of the 
environment, it would then be unusual if z nbi = h(r) . Rather, 
the bottom lies just below the index 




(3.14) 



such that z nbi <h(r). To perform the integrals needed, how- 
ever, the conditions 

zi = 0 , (3.15) 

ZNs = h(r) , 

must be satisfied. This may be accomplished by first trans- 
forming to a zj space (which is not equally spaced) , which 
satisfies Eq. (3.15), and then interpolating Zj to z^ where 



Zk 




(3.16) 



Part of the difficulty is to define the sound speed on 
the non-grid points z = 0 and z = h(r) . Fortunately, this is 
essentially already accomplished in the MMPE model. The 
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surface sound speed is required as an input to the MMPE, and 
is stored in the sound speed array by the environmental prop- 
agator subroutine. The sound speed at the bottom depth, 
h(r) , is also already estimated at the grid point nbi . While 
this is an estimate, it provides a good approximation for 
cases where energy doesn't significantly interact with the 
bottom. This results in 

Z\ = 0=> z = 0 => c(0) = SS , (3.17) 

z Ns = h(r)=*z = h(r) =>c(h(r)) = SSBW , 

where SSO is the surface sound speed input into MMPE and SSBW 

is the estimated water sound speed at the water/bottom 

interface . 

The numerical scheme for computing c 0 and the tilde 
transformation can now be defined. Using Eqs . (3.6) and 

(3.12) , 



c 0 (r) = (3.18) 



h 2 (r) 


z i 


1 


" nbi 

y z k _z k-l 


I 


h ( r )-z„bi 1 


4 


J c] + Jsso_ 


i 


_k = 2 ^ + V C k-l_ 


t 


[7ssBW + 7c nbi J 
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where it is assumed that the sound speed has already been 
interpolated onto the gridded mesh z k at values c k . Simi- 
larly, Eq. (3.4) with Z\ = 0 becomes 



Z2=2 JT 0 , (3.19) 

zj = Zj-i+2^ ^ -1== < j= 3...nbi + 1 , 

J c k + J c k-l 

ZNs = ^Ns-1 ~ 

Vssbw + 7^: 

Note that Ns = nbi + 2 since the points z = 0 and z = h have 
been added. As a result of this transformation and the above 
definitions , 



c, = SS, (3.20) 

c Ns = SSBW, 
c j * I, Ns = c j - 1 

Given these relationships, it is now necessary to 
interpolate in depth the tilde domain sound speed onto the 
standard grid. Defining 



Zk = 
k = 






Az , 






it is desired that 



c(z k ) = c(z k ) . 



(3.21) 



(3.22) 
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A simple linear interpolation from c(zj) to c(z k ) could be per- 
formed. This, however, would not preserve the linear gradi- 
ent nature of c(z k ) , and the relation defined by Eq. (3.22) 
would not be satisfied. Instead, assume 

Zj_!<z k <Zj, (3.23) 

and then define the fractional depth as 

Zv Zj _ i 

f k = - - J . (3.24) 

Zj - Zj _ l 

Tappert (1995) shows that the inversion of the mapping 
needed produces 

c k = Cj_ 1 +g k (Cj-Cj_i), (3.25) 

where 

gk = + • (3-26) 

V c j + v c j - 1 



C. ACOUSTIC FIELD GENERATION 



In the MMPE model the source amplitude is defined rela- 
tive to a small finite distance from the source. This is 
done to avoid the singularity at range r = 0 in 



= P 0 ^2v Fe ik 0 r 



(3.27) 
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Therefore, the source is defined such that p(r = R 0 ) = P 0 . The 

reference range is defined to be one meter to remain consis- 
tent with most sonar equations. Using this convention, the 
source level is related to P 0 by 

SL = 201og^j dB re P r R 0 (3.28) 

where P r = 1 pPa at the reference range of R 0 . 

The form of the source field '{'(r = 0, z) must now be 
determined. Rewriting Eq. (3.27) as 

^(r, z)e ik ° r = /l p ( r> z) , (3.29) 

M)A/ K 0 

then 

^(r = 0, z) = lim /jf-p(r, z) , (3.30) 

r — >0 1^0 

where in the vicinity of the source the pressure filed takes 

p 

the form of the spherical Green's function. If |p| = — is 
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required at R = R 0 , 



then let a = P 0 R 0 



P = 



_L_ p ik 0 R 

4nR e 



and therefore 



(3.31) 



where R = Jr 2 + z 2 . The source is represented as a point 
source at (0,z s ) by 
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¥( r = 0, z) = a§(z - z s ) . 



(3.32) 



Integrating Eq. (3.30) over all depths of the real ocean 
results in 



a lim I— 



1 



0 p 0 aJR 0 J 471R 



e ik ° r dz . 



(3.33) 



0 

In the far-field, r»z, R may be approximated such that 

R = Jr 2 + z 2 ~r + ^~, z = z - z„ (3.34) 

2r s 

thus reducing Eq. (3.33) to 




0 



and leaving 



(3.36) 



Performing a Fourier transform of Eq. (3.32) with the 
influence of the image source included, and specifying the 
source in the k 2 -domain yields 



*F(r = 0, k) = a J [5(z- z s ) - 5(z + z s )]e lkzZ dz = -2iasin(k z z s ) . (3.37) 
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Equation (3.37) states that the starting field, in the wave- 
number representation, has constant amplitude modulated by a 
phase due to the interaction of the source and its image. 
This representation is consistent with an omnidirectional 
source that supplies equal energy to all wavenumbers. 

Rather than setting the amplitude of this function to 
unity to ensure equal population of all directions, a smooth 
taper is used at high absolute wavenumber values to limit the 
angular dimension of the source and to reduce sidelobe 
influences. This is necessary due to the wide-angle approx- 
imation being valid to only approximately 40° , and due to the 
restriction on how large the finite FFT will allow k z /k 0 . 

The wide-angle source is modified to produce an accurate 
solution in the far-field (Thomson and Bohun, 1988) by 

( k 2V 1/4 

F(kz) = [ 1- k|| >( W <k o ) - (3 - 38) 



The limits on k 2 are chosen as such since for |k z |>k 0 the 
angles of propagation are imaginary. Therefore, the source 



function is tapered within the limits 



K 

k o 



< 1 . 
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Finally, to account for the half-mesh symmetry of the 



Az 
ik — 

depth grid a phase term in the wavenumber domain of e 2 is 
added. The wavenumber domain starting field for the wide- 
angle point source can then be written as 



¥(r = 0,k z ) = -2i 



-sin(k z ) 



k|V l/4 ik^ z 
e 2 ' 
H) 



(3.39) 



Following Tappert ' s procedure, for the COIPE it is eas- 
iest to generate the starting field in the tilde domain. The 
starting depth, z s , must then be transformed to tilde space. 
Defining 




the source depth in the tilde-domain is then 




where 



(3.40) 



(3.41) 



Uj = f Cj + ( 1 — f) Cj _ , . (3.42) 

Note that Eq. (3.42) performs a linear interpolation of c at 
the source depth. Since, from Eq. (3.22) 
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zj-zj-i = 2 jr 0 -=± 



z j Zj _ i 






(3.43) 



then Eq. (3.41) becomes 

2 V^(Zj- z j- 1 ) f 



z s = Zi -1 + - 



2 7 ^( z s- z j-i) 



+ V c j- 1 ^ + «/ c j- 1 



(3.44) 



z s = Zj _ i + 



2 ^( Z S~Zj-1 ) 



(3.45) 



which is merely Eq. (3.19) rewritten using c s = c(z s ) . Also 
from the expression for the tilde domain transformation. 



= JVn(z', r)d z’ , 



(3.46) 



it follows that 



Therefore, 



"T"Z = Jr\(z~T) . 
az 



— = (n(z,r))- |/2 
dz 



(3.47) 



(3.48) 



J( n (z, r))' 



1/2 dz’ . 
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(3.49) 



Since c(z) = c(z) , then n(z) = n(z) , such that 



z z 

z= f(n(z', r))“ 1/2 dz’ = —J= f7c(z')dz' . (3.50) 

J V c o( r ) J 

o o 

Equation (3.50) could be solved numerically, but it does not 
have the analytical expression the previous transformation 
did. Rather, the previous transformation for source depth 
can be applied to receiver depths at the original gridding of 
the environment. In other words, consider 



z k 




(3.51) 



which is essentially what was obtained in Eq. (3.19), where 
now the tilde depths are those found in the original trans- 
formation vector, Zj for j = 2...Ns-l . Below zn s _i , z^ns-i is 
used. 

To obtain the field, however, it must be noted that 

^(z, r) = [n(z, r)]- 1/4v F(z, r) (3.52) 

or, since n(z) = n(z) , 

[n(z, r)] 4 4 / (z, r) = ^(r) . (3.53) 

In either case, H'(Zj) must first be interpolated from 'F(Zk) . 



Since there is no simple gradient structure for 'F(z) / a 
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numerical scheme must be used. For this case, a typical 5- 
point interpolation scheme is used to obtain 'F(Zj) . If Eq. 
(3.52) is used, then n(z k , r) is evaluated and the values of 

^(Zk,r) = [n(z k ,r )] 4 'F(z j ) (3.54) 

are computed where k = j-1 and c(z k ) is evaluated from the 

original profile. Note that [n(zk, r)] 4l F(Zk) could have been 
computed and then interpolated to Zj . However, this would 
only be an approximate interpolation of n . Therefore, the 
former method is an exact interpolation of n and only 

approximate in 'F , which is expected to be the more accurate 
method. 

As was briefly discussed in the previous chapter, the 
COIPE approach is expected to remove the reference sound 
speed ambiguity present in the WAPE. Removing this sensi- 
tivity should provide more accurate results without the need 
to determine an appropriate reference sound speed for the 
specific problem. In the following chapter the COIPE is 
benchmarked against the SPE, WAPE, and normal mode theory. 
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IV. NUMERICAL RESULTS 



The previous chapters present an outline of the histor- 
ical and theoretical background of the parabolic equation 
approximation to the acoustic wave equation as well as the 
COIPE implementation. In this chapter, the focus shifts to 
the analysis of attempts to improve the accuracy and valid- 
ity of the WAPE approximation. First, the effects of improv- 
ing the numerical accuracy of the WAPE, as implemented by the 
MMPE model, is evaluated. Second, the COIPE implementation 
is evaluated under a number of test cases. Finally, the SPE, 
MMPE, and COIPE are subjected to the normal mode decomposi- 
tion technique previously discussed. 

A. DOUBLE-PRECISION IMPLEMENTATION OF THE MMPE 

While the MMPE has been used successfully under a wide 
variety of test cases and environments for a number of years, 
it was recently tested extensively on a number of specific 
problems. These tests were conducted as part of the Shallow 
Water Acoustic Modeling (SWAM) Workshop held at the Naval 
Postgraduate School in 1999. The SWAM workshop provided a 
set of detailed environmental cases that served as the basis 
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for comparison between the numerous modeling efforts being 
undertaken in the United States and internationally. As 
part of this workshop, the convergence sensitivity of the 
MMPE to the choice of the range step, Ar , was analyzed 
(Smith, 2000) . As a result of this analysis, which is dis- 
cussed in more detail later in this section, an effort to 
develop a more numerically precise MMPE implementation was 
undertaken. Increasing the numerical precision of the model 
would also help answer whether truncation of the field, as 
postulated in the introduction, is a source of error in the 
model . 

In an effort to avoid introducing new errors or anoma- 
lies into the MMPE, no algorithm changes were made to the 
program code itself. Rather, the variables and array decla- 
rations, as well as the corresponding post processing rou- 
tines, were recast as double-precision. The original MMPE 
code was, with the exception of the FFT subroutine, predomi- 
nantly single-precision. One of the concerns raised by this 
procedure was that computational times would increase dra- 
matically. This, however, was not the case. Since the bulk 
of the computational work is done by the FFT subroutine, 
which was already in double-precision format, computational 
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times only increased by three to five percent. These 
increases are deemed relatively insignificant and do not in 
and of themselves hinder the use of the double-precision 
MMPE (DP-MMPE) . 

As a method of comparison with the SPE and existing 
MMPE, the DP-MMPE was tested using one of the SWAM test 
cases. The test case chosen was referred to as SWAM Flat-a. 
The environment has an isospeed water column, c 0 = 1500 m/s 

p 0 = LOO g/cm 3 , overlying a flat bottom. The bottom proper- 
ties are defined every 2 km and are linearly interpolated out 
to a range of 20 km. Compressional attenuation is fixed at 
0.1 dB/X and there is no shear in the bottom. The source is 
at a depth of 30 m and a frequency of 250 Hz CW. The envi- 
ronment is shown below in Fig. 4.1. Fig. 4.2 shows the full, 
CW field at 250 Hz for the source at 30 m as obtained from 
using the MMPE model solution. 

During the theoretical development of the wide-angle PE 
approximation it was noted that the range step size Ar , as 
based on the centered-step scheme, is roughly third order 
accurate. The implications of this fact are that as the 
range step size is decreased, the solution accuracy should 
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improve until a stable solution converges. In fact, most PE 
models do follow this behavior. The PE/SSF does not. 






Figure 4.1. The Flat-a test case environment. 



To illustrate the relation between the range-step size 
and the covergence of the solution, a number of solutions 
were generated for the Flat-a case discussed above. The only- 
parameter varied in any of the cases was the size of the 
range-step. In each of the cases the transmission loss (TL) 
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Figure 4.2. MMPE field solution for Flat-a case. 



for a receiver at 35 m is shown in Fig. 4.3. The upper and 
lower panels of Figs. 4.3 and 4.4 show the first and last 5 
km of the solution, respectively. Note that the solution 
appears to be converging as Ar~5 m. As can be seen in Fig. 
4.4, however, the solution diverges for further decrease of 
Ar ; the solution degrades for values of Ar less than about 
5 m. The solution appears to reach convergence at Ar-A,. 
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4.3. Convergence testing for various range-steps. 
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Figure 4.4. Stability testing for various range-steps. 
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A possible cause of this divergence was postulated to 
be truncation of the field and a lack of numerical precision 
in the declaration of variables in the MMPE computer code. 
In an effort to test this theory the DP-MMPE model was run 
with identical model inputs to those used in the convergence 
and stability testing shown above for the MMPE. The results 
of this testing are shown below in Figs. 4.5 and 4.6. 





Figure 4.5. DP-MMPE convergence testing for various range- 

step sizes. 
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Figure 4.6. DP-MMPE stability testing for various range-step 

sizes . 

Clearly, there does not appear to be any discernible 
improvement in the convergence of the solution through the 
use of the DP-MMPE. It exhibits the same general behavior as 
that of the original MMPE model. A comparison between the 
MMPE and DP-MMPE at the 5 m and 2 m range-step sizes is shown 
below in Fig. 4.7. 
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Figure 4.7. DP-MMPE and MMPE comparison for range-step sizes 

of 5 m and 2 m. 

Referring to Fig. 4.7, it can be seen that for a range 
step size of 5 m the solution produced by the two models are 
indeed overlapping. For a range-step size of 2 m the two 
solutions differ slightly with range, with the DP-MMPE 
appearing to diverge less at some range points. One may con- 
clude, however, that while the DP-MMPE does no worse, it does 
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not offer a solution as to why the PE/SSF algorithm seems to 
have a lower limit on the size of the range- step. 

As a result of separate analysis (Smith, 2000), the 
primary cause of this non- convergence was found to be 
related to the structure of the propagator functions. 

Recall that the environmental propagator function is of the 

form e -lk ° ArU °p( r ’ z) / an £ the wavenumber propagator function is of 
the form e - lk o ArT o P ( k z) . The wavenumber propagator decays expo- 
nentially beyond k z >k 0 . From this analysis it can be shown 

that the PE/SSF algorithm, for large range-step size, 
attempts to include a large amount of phase information into 
each range-step. If Ar is chosen to be too large, then 
errors are generated in the solution. Conversely, as Ar is 
reduced in size there is inadequate phase information con- 
tained in each range-step and the solution once again 
becomes inaccurate. The most accurate solutions may be 
expected for range-steps where a full cycle of phase infor- 
mation is contained in each propagator. The level at which 
Ar appears to provide the greatest degree of convergence is 
Ar~X (Smith 2000) for shallow water environments such as 
tested here. 
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Before leaving the area of solution convergence based 
on the range- step size, a final comparison between the MMPE 
and the COIPE is made. While the COIPE model was not imple- 
mented specifically to address this problem, it is of inter- 
est as to whether any of the changes made in the COIPE affect 
the convergence due to Ar . Below, in Fig. 4.8, the COIPE 
model is compared against the MMPE for the same test case 
used above. The MMPE results were chosen over the DP-MMPE 
since the DP-MMPE solutions did not offer any definitive 
advantages over the normal implementation. As Fig. 4.8 
shows, the COIPE does not appear to offer any greater degree 
of convergence for the given Ar sizes in this environment. 
These results are inconclusive in regards to the COIPE 
model, but the implementation was not expected to affect the 
Ar sensitivity. 

Finally, as a means of comparison against another mod- 
eling technique for this test case, the MMPE results for 
range- step size 5 m are compared to those developed by Mikhin 
(2000) for SWAM. These results were chosen as a benchmark 
due to the high degree of agreement with several test cases 
and several modeling techniques. 
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Figure 4.8. COIPE and MMPE convergence testing. 

The comparison of the MMPE solution with the solution 
provided by Mikhin is shown below in Fig. 4.9. Solutions for 
rang-step sizes of 20, 5, and 2 m were chosen to demonstrate 
the degree of convergence with the reference solution for 
varying Ar . The figure shows very good agreement between 
the MMPE and Mikhin solutions at short ranges. The two tech- 
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niques show a greater degree of separation with increasing 
range. This is most likely the result of cumulative phase 
errors inherent to the wide-angle PE technique employed 
here. The figure also demonstrates, once again, that the 





Figure 4.9. MMPE and Mikhin solution comparison for various 
range step-sizes. 



choice of a 5 m range-step size provide the greatest degree 



of agreement with the reference solution. 
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B. THE COIPE IMPLEMENTATION OF THE MMPE 



As was discussed in the last chapter, the COIPE model 
was developed by Tapper t (1995) in an attempt to remove or 
reduce the dependence of the wide-angle PE approximation on 
the choice of the reference sound speed value while increas- 
ing the overall accuracy of the model. While this sensitiv- 
ity has been observed in a number of test cases, it was 
extensively studied during one of the PE Workshop II (PE- II) 
test cases (Porter and Jensen, 1993) . This case, referred to 
as the Porter duct problem, was based on the observation that 
many of the wide-angle split- step PE approximations overes- 
timated the value of the TL when applied to long-range prop- 
agation in a leaky surface duct environment. The SPE 
implementation of the PE/SSF algorithm does not generate 
this anomaly. Additionally, other PE implementations based 
on Pade' series expansions of the square root operator 
(which maintain the proper normal mode basis set) do not 
result in this overestimation of the downrange TL. The 
results of PE-II showed that many wide-angle models that do 
not use the Thoms on -Chapman (TC) wide-angle approximation 
avoid this problem. 
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The Porter duct environment is one of a source and 



receiver located in a 250 m deep surface duct. The source is 
at a depth of 25 m and frequency of 80 Hz. The receiver is 
at a depth of 100 m. The ocean bottom is range- independent 
with a lossy fluid half -space beginning at a depth of 4 km. 

Of the 78 propagating modes in the water column, one is 
trapped in the surface duct. Figures 4.10 and 4.11 illus- 
trate the sound speed profile and TL field for this environ- 
ment respectively. 




Figure 4.10. Porter duct sound speed profile plot. 
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Figure 4.11. TL field plot for the Porter duct problem. 

Finn Jensen, during PE Workshop II, showed that the SPE 
was in excellent agreement with the normal mode reference 
solution used in the test case (Porter and Jensen, 1993). 
Because this reference solution is no longer available, the 
SPE results were used as the reference for this analysis. 
The only alterations to the various model input parameters 
changed in the following analysis is that of the reference 
sound speed, c 0 . Prior to looking at the COIPE model or the 
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effects of varying c 0 , the TL plot for the receiver depth of 
100 m is shown in Fig. 4.12 for the SPE and WAPE implementa- 
tions when a reference sound speed of c 0 = 1500 m/s is used. 

It can be clearly seen that the WAPE solution shows a dis- 
tinct drop off at ranges beyond about 60 km. It is this 
anomaly that the COIPE model is meant to eliminate. 




Figure 4.12. TL Plots for SPE and WAPE at the receiver depth 
of 100 m and c 0 = 1500 m/s. 

The reference sound speed was then set at 1485 m/s and 
the two models were run again. These results are shown in 
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Fig. 4.13. This plot illustrates the effect the value of c 0 
has on the WAPE solution. By "correctly" choosing the value 
of c 0 / the WAPE correctly determines the long range TL. The 
difficulty, however, is in correctly determining this value. 
A change in c 0 by as little as 5 m/s 

to either side of this "correct" value results in the long 
range drop off shown in Fig. 4.12. 




Figure 4.13. TL Plots for SPE and WAPE at the receiver depth 
of 100 m and c 0 = 1485 m/s. 



66 




The COIPE was then tested using identical environmental 



inputs to those used in the analysis above. COIPE model runs 
were made at c 0 values of 1500 and 1485 to provide comparison 
with the SPE results. These are shown in Figs. 4.14 and 4.15 
respectively. The COIPE model compares well with the SPE 
solution at all ranges and there is no drop off at extended 
range as was observed with the WAPE. The COIPE solution is 
indeed insensitive to changes in the reference sound 




Figure 4.14. TL Plots for SPE and COIPE at the receiver depth 
of 100 m and c 0 = 1500 m/s. 
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Figure 4.15. TL Plots for SPE and COIPE at the receiver 
depth of 100 m and c 0 = 1485 m/s. 

speed, as shown in Fig. 4.15. While both plots in Fig. 4.15 
do show phase differences from those of Fig. 4.14, the COIPE 
follows the reference solution equally well in both cases. 
The COIPE was tested at a wide variety of c 0 values with sim- 
ilarly successful results. 

Clearly the COIPE corrects the TL drop-off encountered 
in the MMPE and similar wide-angle implementations. The 
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COIPE corrects the small error in the phase calculation 
caused by the incorrect choice of the reference sound speed. 

A small error in this phase calculation can cause a large, on 
the order of 20 dB, downrange TL error (Porter and Jensen, 
1993) . The test case highlights this problem because it sup- 
ports a single trapped mode in the surface duct that is a 
leaky mode. The trapped mode continually loses energy into 
the lower region where, due to the strong upward refracting 
profile, it is directed back up toward the surface duct caus- 
ing interferences. Small phase errors result in incorrectly 
applying these interferences to the surface duct mode. The 
transformation and subsequent calculation of the reference 
sound speed at each range step removes the source of these 
phase errors and thus correctly applies the interferences 
described above. 

C. NORMAL MODE DECOMPOSITION OF THE PE FIELD 

The mode functions of the wide-angle PE approximation 
form a different basis set for modal expansion than those 
obtained using standard normal mode theory. As discussed in 
the theoretical development presented earlier, the SPE does 
reproduce the standard basis set. This section presents an 
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analysis and comparison of the modal decomposition of the 
SPE, WAPE, and COIPE models. 

The ocean environment chosen for this analysis is the 
Munk canonical profile (Munk, 1974) . This is a deep ocean, 
range- independent environment with sound speed profile char- 
acterized by 

[c(z)-c 0 ]/c 0 = e(e^ -r| - 1) , (4.1) 

where 8 = 0.0057, the scaled depth variable T) = 2(z - z axis )/B , 
and the reference sound speed is 1490 m/s. The axis depth 
was set at z axis = 1000m, and the bottom depth at 4500 m. The 
sound speed profile is shown below in Fig. 4.16. The source 
was placed at a depth of 1000 m with a frequency of 100 Hz. 

A computational range of 100 km was used for PE field calcu- 
lations. A bottom depth of 4500 m and a computational depth 
of 8000 m were used. 

The acoustic pressure field was calculated using each 
of the three models. The field is shown in Fig. 4.17. The 
resulting fields were then decomposed into normal modes and 
the mode amplitudes determined using the techniques outlined 
in Chapter II for each of the models. Since the 
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Figure 4.16. Munk canonical sound speed profile, 
environment chosen is range- independent , the mode amplitudes 
should remain constant with range. The amplitudes for 
selected modes are shown in Figs. 4.18 thru 4.20 for the SPE, 
WAPE, and COIPE respectively. Modes were chosen to avoid 
overlapping plots to improve readability and to avoid bottom 
interactions. The reference sound speed is 1500 m/s for this 
series of plots. Each of the three models produce satisfac- 
tory results for the lower modes. Note the nearly constant 
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Figure 4.17. TL Field plot for Munk canonical profile, 
amplitude with range for all three models. The higher order 
modes of the WAPE plot, however, show increasing levels of 
fluctuations. Conversely, the COIPE modal amplitudes com- 
pare well with the SPE results and do not show the high 
degree of fluctuation seen in the WAPE plots. 

The three models were then tested for this environment 
using a reference sound speed of 1490 m/s. The results are 
shown in Figs. 4.21 through 4.23. The results for the COIPE, 
Fig. 4.23, are visibly unchanged, as expected. The WAPE 
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Figure 4.18. SPE mode amplitudes with range for modes 1, 5, 
10, 15, 20, 25, 30, 35, 40, 45, and 50. (c 0 = 1500 m/s) 




Figure 4.19. WAPE mode amplitudes with range for modes 1, 5, 
10, 15, 20, 25, 30, 35, 40, 45, and 50. (c 0 = 1500 m/s) 
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Figure 4.20. COIPE mode amplitudes with range for modes 1, 5, 
10, 15, 20, 25, 30, 35, 40, 45, and 50. (c 0 = 1500 m/s) 



Figure 4.21. SPE mode amplitudes with range for modes 1, 5, 
10, 15, 20, 25, 30, 35, 40, 45, and 50. (c 0 = 1490 m/s) 
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Figure 4.22. WAPE mode amplitudes with range for modes 1, 5, 
10, 15 , 20, 25, 30, 35, 40, 45, and 50. (c 0 = 1490 m/s) 
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Figure 4.23. COIPE mode amplitudes with range for modes 1, 5, 
10, 15, 20, 25, 30, 35, 40, 45, and 50. (c 0 = 1490 m/s) 
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amplitudes show somewhat less fluctuation for the higher 
modes. This can most likely be attributed to the sensitivity 
to the reference sound speed discussed in the previous sec- 
tion. 

This analysis has demonstrated the wide-angle PE's lack 
of the ability to reproduce the standard normal mode basis 
set as defined by normal mode theory and as duplicated by the 
SPE. The WAPE shows increasing levels of fluctuation in the 
modal amplitudes with increasing mode numbers. This fluctu- 
ation results in the produced PE field not exactly repre- 
senting the environment under study. The COIPE, however, 
does not show the same degree of fluctuation in modal ampli- 
tudes, and therefore should better represent the acoustic 
field. It should be noted that the modal amplitudes of the 
COIPE do not exactly match those of the SPE. This is due to 
the generation of the source function in the tilde-domain. A 
more proper treatment might generate the starting field in 
real space and then transform to tilde space. However, this 
would introduce some numerical errors to the transformation 
of the filed form tilde to real space. This is an area for 
further analysis and study. 
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V. CONCLUSIONS AND SUMMARY 



This thesis has examined the parabolic equation approx- 
imation to the acoustic wave equation and, in particular, 
two possible implementations aimed at removing inherent WAPE 
errors were analyzed. The theoretical background underlying 
each of these models was presented along with the specific 
numerical implementation used for testing and comparison 
with existing PE approximations. Finally, numerical results 
were provided comparing the SPE, WAPE, DP-WAPE, and COIPE 
models against a number of test cases. These included 
results from the SWAM workshop, the Porter duct problem, and 
the Munk canonical profile. 

The parabolic equation and many of its derivatives suf- 
fer from a number of anomalies or errors depending on the 
particular implementation being used. Of these, three were 
examined in this analysis. First, the selection of the 
range-step size used by the PE/SSF algorithm affects both 
the convergence and the stability (Smith, 2000) of the solu- 
tion with an error to second or third order in the range-step 
increment (Jensen, et . al., 1994). Second, the WAPE suffers 
an anomaly that under certain environmental conditions it 
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incorrectly computes the down- range TL (Porter and Jensen, 
1993) . This anomaly is manifest in a distinct drop in TL. 
Results of the Parabolic Equation Workshop II (PE- II) 
attribute this error to incorrectly handling the loss of 
energy from a leaky surface duct mode and its subsequent 
interaction with the remaining sound field. It was shown 
here, and in PE- II, that the manifestation of this anomaly is 
highly dependent on the choice of the reference sound speed. 

A small change in the value of c 0 results in a significant 
change in the computed down-range TL. Third, the WAPE does 
not reproduce the standard normal mode basis set as defined 
by normal mode theory and as reproduced by the SPE. In addi- 
tion to the phase errors experienced to some degree by all PE 
approximations, the WAPE mode amplitudes exhibit a range 
dependence in a range- independent environment (Smith and 
Smith, 1997). This effect was found to contribute to the 
anomalous TL findings just described. 

A double-precision version of the WAPE model was imple- 
mented in an attempt to remove or improve the convergence 
sensitivity to the selection of the range-step size used by 
the PE/SSF algorithm. The double-precision model was tested 
against the existing SPE, WAPE, and COIPE models in a simple 
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range-dependent environment originally defined for SWAM. 

The double-precision model was found to do as well as, but no 
better than, the existing approximations under the tested 
conditions. From this analysis it was concluded that the 
range-step convergence sensitivity is inherent to the PE/SSF 
algorithm used by the models. This is supported by separate 
analysis done by Smith (2000) . The most accurate solutions 
may be expected for range- steps where a full cycle of phase 
information is contained in each of the propagators used by 
the algorithm. The level at which Ar appears to provide the 
greatest degree of convergence is Ar-X (Smith 2000) for 
shallow water environments such as those tested. The COIPE 
model was also tested for sensitivity to the choice of the 
range-step on solution convergence. It was found to be sim- 
ilarly affected, as was expected. 

The COIPE model, a WAPE implementation, was developed 

to remove or reduce the solution sensitivity to the choice of 

the reference sound speed value (Tappert, 1991) . The COIPE 

model uses an operator transformation which results in a 

range -dependent series of values for the reference sound 

speed. The COIPE model was compared with the SPE and WAPE 

models using the Porter duct problem developed for PE- II 
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(Porter and Jensen, 1993) . These tests were conducted using 
a variety of reference sound speed values, with successful 
results in all cases. The COIPE model does not exhibit the 
TL drop with range seen with the WAPE. 

Finally, the COIPE was tested against the SPE and WAPE 
models for its ability to decompose properly into the stan- 
dard mode basis set. The models were tested in a deep-ocean 
Munk canonical profile environment. As discussed, the WAPE 
does not correctly decompose with this basis set, and shows a 
range dependence in its modal amplitudes in a range- indepen- 
dent environment. This range dependence can be seen as fluc- 
tuations in the modal amplitudes with range. The COIPE modal 
amplitudes show much less fluctuation with range than those 
of the WAPE, and more closely reproduce those obtained by the 
SPE. This is a result of the improved handling of the phase 
by the COIPE as outlined above and demonstrated through the 
modal decomposition. 

The double-precision model testing provides further 
evidence, in addition to that shown by Smith (2 000) , that the 
small range-step sensitivity is a fundamental property of 
the PE/SSF algorithm, and no further investigation of this 
property is necessary. The results of the COIPE model test- 
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ing, while promising, are not conclusive. The COIPE cor- 
rects the error seen in the Porter duct problem, and improves 
the reproduction of the standard normal mode basis set, but 
requires further investigation. Preliminary testing in 
shallow water environments indicate the COIPE may not pro- 
duce results significantly different from those of the WAPE 
model. Additionally, the mode amplitudes themselves, while 
showing less fluctuations, are of different relative magni- 
tudes from those produced by the SPE. This remains as an 
area for further investigation. 
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